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I. INTRODUCTION 



The purpose of this paper is to develop analytic expressions for option-pricing of markets 
with fluctuating volatilities of a given distribution. There are several good reasons for con- 
sidering the statistical properties of volatilities. For instance, a number of important models 
of price changes include explicitly their time-dependence, for example the Heston model [lj], 
or the famous ARCH j^, GARCH 0], and multiscale GARCH [4] models. Changes in the 
daily volatility are qualitatively well explained by models relating volatility to the amount 
of information arriving in the market at a given time 0] . There is also considerable practical 
interest in volatility distributions since they provide traders with an essential quantitative 
information on the riskiness of an asset [H, E}. As such it provides a key input in portfolio 
construction. 

Unlike returns which are correlated only on very short time scales [sj of a few minutes and 
can roughly be approximated by Markovian process, the volatility changes exhibit memory 
with time correlations up to many years U 0, Q E3- In Ref. [ll| it was shown that 



the Standard & Poor 500 (hereafter S&P 500) volatility data can be fitted quite well by a 
log-normal distribution. This result appears to be at odds with the fact that the log-normal 



shape of the distribution typically signalizes a multiplicative nature [13j of an underlying 
stochastic process. This is rather surprising in view of efficient market theories [si which 
assume that the price changes are caused by incoming new information about an asset. Such 
information-induced price changes are additive and should not give rise to multiplicative 
process. We cure this contradiction by observing that the same volatility data can be fitted 



equally well by a Chi distribution [141 . Il5l |. In Appendix A we show that corresponding 
sample paths follow the additive rather than multiplicative Ito stochastic process. With the 
help of Ito's lemma one may show (cf. Appendix A) that the variance v it) = a 2 (t) follows 
the Ito's stochastic equation 



dv(t) = 7(*)[f(*) ~M*M*) -a(v(t),n(t),v(t))]dt + y/2^(t)v(t) dW(t) . (1) 



where W(t) is a Wiener process. Here , y(t),n(t) and u(t) are arbitrary non-singular positive 
real functions on Function a(. . .) is non-singular in all its arguments and it tends to 
zero at large £'s. From the corresponding Fokker-Planck equation one may determine the 
time-dependent distribution of v which reads 

Ut)Ht){v) = j^ry W)\ V(t) v with J dv f^v) = 1. (2) 

The distribution fn,v(v) is the normalized Gamma probability density [14j, whose profile 
is shown in Fig. El It has an average v = v/fi, a variance (v — v) 2 = v/fi 2 , a skewness 
(v — v) 3 = 2/ yjv, and an excess kurtosis (v — v) 4 / (v — v) 2 — 3 = 6/V (see, e.g., Refs. (l4l. [Hi] ) . 
The Gamma distribution (|2J) will play a key role in the following reasonings. The Chi 
distribution p(a, t) is related with the Gamma distribution through the relation 

p{a,t) = 2af fl{t)Mt) {a 2 ) . (3) 
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Often is the functional form ([3]) itself called a Gamma distribution. To avoid potential 
ambiguities, we shall confine in the following to the name Chi distribution. The derivation 
of p(a, t) from the underlying additive process, rather than a multiplicative process, shows 
that the Chi distribution is well compatible with efficient markets. 

The Gamma distributions for fluctuations of v allow us to generate an entirely new 
class of option pricing formulas. For this we use the well-known fact of non-equilibrium 
statistical physics 171 ]. that the density matrix of a system with fluctuating temperature can 
be written as a density matrix for the system with fixed temperature averaged with res pec t 



to a temperature distribution function. Path integrals conveniently facilitate this task 16 



With the help of the so-called Schwinger trick we show that if the distribution of the inverse- 
temperatures of the log-returns is of the Gamma type, the distribution in momentum space 
is of the Tsallis type. To put this observation into a relevant context, we recall that Tsallis 
distributions in momentum space enjoy a key role in statistical physics as being optimal 
in an information theoretical sense: given prior information only on the covariance matrix 
and a so-called escort parameter, they contain the least possible assumptions, i.e., they are 
the most likely unbiased representation of the provided data. Some background material on 
this subject is reviewed in Appendix C. The observed connection with Tsallis distribution in 
momentum space can be fruitfully used to address the issue of the time-compounded density 
function for stock fluctuations with Chi distributed volatility. The Markovian property of the 
market prices then simply translates to product of distributions in momentum space. At the 
level of the distributions of log-returns, this automatically gives rise to the correct Chapman- 
Kolmogorov equation for Markovian processes. The time-compounded distribution function 
constructed by convolution of integrals represents the desired measure of stock fluctuations 
at given time t. 

It should be emphasized that since the distribution of log-returns is related to the Tsallis 
distribution in momentum space by a Fourier transform, the log-return data do not inherit 
the heavy tails of the Tsallis distribution in momentum space. In this respect our approach 
is quite different from data analysis based on a Tsallis distribution of log-returns, such as 



the so-called 'non-extensive thermostatic" 18] . The formal advantage of our model lies in 



its optimal use of the mathematical machinery of the Black-Scholes theory. In particular, 
one still has the option pricing formula that is linear in the spot probability of the strike- 
price payment. Put and call options still obey the put-call parity relation, and the A- hedge 
still coincides with the spot probability that multiplies the asset price in the Black-Scholes 
formula. On the other hand, the desired features such as peaked middles of financial asset 
fluctuations and semi-fat tails are present for sufficiently long times. The latter ensures that 
ensuing option prices can differ noticeably from the Black-Scholes curves for rather long 
expiration dates (days or even moths) despite the validity of central limit theorem (CLT). 
One may thus expect our formulation to be useful, e.g., for short or mid-term maturity 
options ( "mesoscopic" time lags). 

The paper is organized as follows. In Section [III we fit the high-frequency data set of the 
volatility fluctuations in the S&P 500 index from 1985 to 2007 to a Chi distribution. We 
also show that this fit holds well for different time averaging windows. In Section II III we 
represent the distribution of log-returns as a superposition of Gaussian distributions with 
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the above volatility behavior and find that the associated distribution in momentum space is 
of the Tsallis type. By assuming further that the compounded stock fluctuations over longer 
times are Markovian we obtain in Section HVl the corresponding natural martingale measure. 
With it we derive an option-pricing formula as a superposition of Black-Scholes expressions. 
The departure from Black-Scholes results is expressed as a sum of two qualitatively differ- 
ent expansions: expansion in powers of the excess curtosis and expansion in momenta of 
moneyness. The former represents the expansion in market characteristics while the latter 
is an expansion based on contract characteristics. In Section [V] we determine the crossover 
time below which our option prices differs from Black-Scholes. Comparisons with observed 
log-returns are presented in Section IVII There we show that for the Dow Jones Euro Stoxx 
50 data the amount of the residual risk in the A-hedge portfolio is very small and that the 
crossover time is roughly 7 — 8 months. Section IVIH finally, is devoted to conclusions. For 
reader's convenience we also include three appendices. In Appendix A we solve and discuss 
the Fokker-Planck equations that are associated with stochastic equations for the volatility 
and variance presented in Introduction. In Appendix B we deal with some mathematical 
manipulations needed in Section IIVI In Appendix C we present some basics for Renyi and 
THC statistics to better understand some remarks in the main body of the paper. 



II. EMPIRICAL MOTIVATION 

Our work is motivated by data sets on volatility distributions extracted from the S&P 
500 stock market index. We briefly remind the reader the procedure for obtaining these 



numbers. A detailed description can be found in Ref . [1 1| and citations therein. 

It is well known that the autocorrelation function of stock market returns decays expo- 
nentially with a short characteristic time - typically few trading minutes (e.g., ~ 4 min. for 
S&P 500 index). Hence the log- returns R are basically uncorrelated random variables, nev- 
ertheless, they are not independent since higher-order correlations reveal a richer structure. 



In fact, empirical analysis of financial data confirms [llj that the autocorrelation function of 
non- linear functions of R, such as \R\ or R 2 , has much longer decorrelation time (memory) 
spanning up to several years (few months for S&P 500 index). These observations imply 
that a realistic model for the return-generating processes should account for a non-linear 
dependence in the returns. The latter can be most naturally achieved via some additional 
memory-bearing stochastic process. Most simply, this subsidiary process involves directly 
the volatility. 

One may bypass a construction of the volatility model by trying to define a "judicious" 
volatility estimator directly from financial time series. Difficulty resides, however, in the 
fact that although the volatility is supposed to be a measure of the magnitude of market 
fluctuations, it is not immediately clear how such a measure should be quantified. Among 
many definitions present in the literature [l9] we focuss here on the estimator proposed in 
Ref. [1 1| . There one defines the volatility ay at a given time t v as the arithmetic average 
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of the absolute value of the log-returns 



R{t r , 



hi 



S{t n + At) „ S{t n + At)-S{t n ) 



S(t n ) S(t n ) 
over some time window of the length T = nAt (n G N + , At is the sampling interval), i.e. 



(4) 



n— I+J7 



(t^ = r/ At, t m = mAt) . 



(5) 



m=r\ 



Because (jSJ) is basically a forward-time mean value of absolute returns, it may indeed serve 
as a good measure of market fluctuations. 

With this 0r(t) we can construct the corresponding volatility probability density function 
according to a relative frequency prescription: 



Pt(c")A(T 



# a T (t) E [a,a + Aa] 



n 



(6) 



In the limit of the very large time window T, one might expect Pr (o') to approach a Gaussian 
distribution, since the CLT holds also for correlated time series [20(], although with an often 
slower convergence than for independent processes 
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There is, however, a logical caveat in the definition (jSJ). The observed volatility fluctua- 
tions have a spuriously superposed pattern caused by intra-day fluctuations. This is because 
over the day, the market activity is large at the beginning and at the end, but exhibits a 
broad minimum around a noon 22] . Since the volatility is supposed to measure the magni- 
tude of the market activity, this intra-day pattern should be removed in order to avoid false 



correlations. This can be remedied by considering the normalized returns 11 



\n[S(t n + At)/S(t n )} „ S(t n + At) - S(t n ) 



(\\n[S(t n + At)/S(t r 



S(t n ) 



S(t n + At) - S(t r 



S(tr. 



■ (7) 



Here (. . .) is the average over the whole sequence of trading days, i.e., 



, iv td 

(\\n[S(t n + At)/S(t n )}\) = —Y^lHSjitn + A^/Sjit 



(8) 



with Sj(t n ) denoting the spot price of the index at the time t n on the jth trading day, and 
N td denotes the total number of trading days. So in (0) each spot return is divided by its 
natural "spot return scale" . The corresponding normalized volatility is then defined as 



n— 1+7) 



C"Tnor(^r)) = — / , \Rnor{t 

n ^-^ 



(9) 



m=r\ 
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This is the quantity whose distribution will be discussed here (we omit the subscript "nor" 
in the following). To this end we shall examine the data set for the S&P 500 stock market 
index. Data in question were gathered over the period of 22 years from Jan 1985 to Jan 
2007 at roughly 5 minute increments. The corresponding empirical time sequence is seen in 
Fig. [[J With the help of the prescription (jSJ) we obtain the volatility for the above S&P 500 



i ' i • i ' i ' r 




, i , i i i , i , i 

5/1/1985 5/1/1989 5/1/1993 5/1/1997 5/1/2001 5/1/2005 

Time 



FIG. 1: Logarithmic plot of the S&P 500 index S(t) over 22 years (5 January 1985 - 5 January 
2007) with sampling intervals At = 5 min. The linear fit shows the typical exponential growth at 
an annual rate of ~ 15%. Only end-of-day prices are shown. 

index data shown in Fig. [2j The corresponding normalized volatility ar nor (t) is shown in 




FIG. 2: Volatility a T (t) of S&P 500 index over 22 years (5 January 1985 - 5 January 2007) with 
sampling intervals At = 5 min. and time window T = 300 min., i.e., roughly one trading day. 
Only 3000 volatility values are plotted. 



Fig- El The normalization was taken with respect to N t & = 5550 (i.e., approximate number 
of trading days between 5 Jan. 1985 - 5 Jan. 2007). Figure H] shows the volatility probability 
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5/1/1985 5/1/1989 5/1/1993 3/1/1997 5/1/2001 18/1/2005 



Time 



FIG. 3: Normalized volatility o"Tnor(£) of S&P 500 index over 22 years (5 January 1985 - 5 January 
2007) with sampling interval At = 5 min., Ntd = 5550, and time window T = 300 min., i.e., 
roughly one trading day. Only 1000 volatility values are plotted. 



density function pr{c) with the time window T = 300 min. While it is seen that pr{c) can be 

), the Chi distribution gives 



well fitted with the Log- normal function (as proposed in Ref . [11 
a better fit in the central part and is of roughly the same quality in the tail part. Because the 
value of the volatility quantifies the asset risk, Fig. H] implies that the Log- normal fit slightly 
overestimates, while the Chi distribution slightly underestimates the large risks. In addition, 
the functions fi and v that parameterize the Chi (and Gamma) distribution (cf. Eq.([3])) are 
in the observed time windows proportional to T, i.e. fi(T) = \xT and v(T) = vT . This 
tendency is confirmed in Fig. [5] which shows the corresponding variance distribution pr{v) 
for three different window sizes T. For better comparison we use the scaled distribution 




h Empirical volatility 

Gauss distr. fit 

Log-normal distr. fit 

Chi distr. fit 



Normalized volatility 



FIG. 4: Best Gaussian, Log-normal and Chi-distribution fit for the S&P 500 empirical volatility. 
The time window T = 300 min. 



form, T(Tu)pr(v)/(T}i) as a function of T\iv + (1 — Tv) \n(Tfw), with \i = E(v)/Var(v) 
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and v = E 2 (v)/Var(v) (E and Var are empirical mean and variance, respectively). After 
the above scaling, empirical variance distributions, with T = 300, 600 and 900 min. nicely 
"collapse" to the canonical exponential exp(— x) (i.e., Gamma distribution with both mean 
and variance equal to 1), confirming thus the assumed Gamma distribution behavior for 
rather long averaging times. In the limit of very long T, one expects that Pt(c") becomes 



, A 



• T = 900 min. 

■ T = 600 min. 

A T = 300 min. 

— exp(— Vtr) 



•■A 



n& A A A *--A -A- 4BA - A- A«A- A— * -A» A— — A A- 



Transformed variance vtr 



FIG. 5: The scaled variance distribution for different window sizes T. The scaled distribution is 
plotted as a function of the transformed volatility vt r = Tfiv + (1 — Tv) ln(Tfj,v). 



Gaussian, due to CLT. For the times considered here, however, the Chi distribution fits 
the data better than a Gaussian one. We see thus that the long-range correlations in 
the volatility fluctuations considerably slow down the convergence towards the Gaussian 
distribution. The conclusions that one my draw from this is that, even for large times, 
the tail of pr(c") has a fairly large amount of distribution (or information) that cannot be 
ignored. This will provide a key input in generalizing the Black-Scholes pricing formula in 
Section HVl 



III. THEORETICAL DIGRESSION TSALLIS' DENSITY OPERATOR 



We have mentioned before the result in statistical physics [r7[ that a density matrix of 
a system with fluctuating temperature can be written as a density matrix with fixed tem- 
perature averaged with respect to some temperature distribution function. This technique 



was recenly proposed as a natural frame for composing non-Gaussian price distributions 23 
and for deriving option pricing formulas 



24 



Consider the Gaussian distribution of the Brownian motion of a particle of unit mass as 
a function of the inverse temperature (3g = l/T: 

p G (x b ,x a ,p G ) = — L= e-to-" )2/2fe • (10) 
V2vr ( £J G 
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We use natural units in which the Boltzmann constant ks has the value 1. The variance v = 
a 2 of this distribution is obviously equal to (3g- Following Ref. 23] we form a superposition of 



these distributions as an integral over different inverse temperatures, i.e., different variances 
v: 

roc i 

Ps (x b ,x a ,P) = / dv f x/s ( v )-— e -to— /3 = v = u/fi, (11) 
Jo \Jlnv 

where f^ v {y) is the Gamma distribution defined in (j2J) whose average lies at v — vj\i. The 
subscript 5 of the distribution characterizes the ratio 



v v 2 

a quantity which we shall call the spread of the Gamma distribution. In terms of the 
integration variable s = piv, we can rewrite (fTT!) in the form 

rt^fl) = rmfT s " ,e "^ ^ fm ' ■ < 13 ' 

In quantum mechanics, one writes such a distribution in a notation due to Dirac as a matrix 
element 

p(x b ,x a ;p) = (x b \p(p)\x a ) , (14) 



of a density operator 



m = 



1 /'■ x d£ a i/« e -. e -v^ f /9 = 1/ / /4= i/ /l< y ! 



r(i/5) A * 



(15) 



where 



i? = 



P 



2 d 2 



(16) 



is the Hamilton operator of a freeparticle of unit mass. The integral over s can be done (as 
in the so-called Schwinger trick [16j) and yields: 



p((3) = 1 + pSH 



-1/8 



(17) 



This is the (un-normalized) Tsallis density operator (cf. Appendix B) with a so-called 
escort parameter q related to the spread parameter 5 by q = 1 + 5, and an inverse Tsallis 
temperature (3 equal to v j p. In the limit 5^0 where q — ► 1 + , the Tsallis operator (|T7|) 
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FIG. 6: Gamma distribution for various values of ji and u, with average v = u/fi = (3, variance 

2 

(v — v) 2 = vj '/i 2 , skewness (v — v) 3 = 2/y/u, and excess kurtosis (v — u) 4 /(f — v) 2 — 3 = 6/za 



converges to an exponential 

m - e-*> A , Po = (3\ q=1} (18) 



and the distribution function (|T3|) becomes Gaussian. The corresponding limit of the Gamma 
distribution (|2J) is a 5-function: 

Hence the density matrix (11 ip approaches for 5 — > the canonical density matrix of the 
Gibbs-Boltzmann statistics of the inverse temperature v j \i = fla- Note, that in order to 
ensure that (3q is finite in the small-5 limit, [i must behave as 1/50G at 5 — > 0. The 
reader may easily check that a superposition of the type (fTTl) exists also for distributions of 
Bashkirov's 1-st version of thermostatics (cf. Appendix B). In that case 0<q< 1,5 = 1 — q 
and j3 i — {3 . In the limit q — > 1_ the density matrix becomes again a Gibbs-Boltzmann 
canonical density matrix. 



IV. GENERALIZED OPTION PRICING FORMULA 

Consider a continuously tradable stock and assume that the stock fluctuations over short 
time intervals At® such as Ato = 1 day are described by the stationary Tsallis density 
operator (|17p . The time Ato will from now on play the role of a time unit. We now analyze 



the modifications of the Black-Scholes formula for pricing European call options 25] brought 



about by the fluctuations of the volatilities. After a time t > (always in units of At ) the 



10 



stock fluctuations follow a Tsallis density operator 



1 + 35H 



-t/6 



17 -t poo J 

' _ s t/6 e -s[l+f35H] 



t \ t/s z-* 
P5J rjt/s) 



dv 



r (t/5) 



v t/5 e -tvn e -tvH 



(20) 



where Z is a normalization factor. In the sequel we shall allow for a drift of the returns by 
extending the Hamiltonian operator (1161) to 

H = p 2 /2+pr Xw /v, (21) 

where r\y is the growth rate of the riskfree investment, and r Xw = r-w + v/2 is the associated 
growth rate of its logarithm. The parameter j3 equals to 1/Sfi, as before. The matrix 
elements of (12"0"|) yield the time-compounded probability density 



P s {x b ,t b ;x a ,t a ) = (x b |[p(/?)]*|x a ) , t = t b -t a > 



(22) 



Being the matrix element of the product of operators p(/3). Ps(x b ,t b ; x a ,t a ) fulfills trivially 
the Chapman-Kolmogorov relation for a Markovian process 

/oo 
dx P s (x b ,t b ;x,t c )Ps(x,t c ;x a ,t a ) , t b > t c > t a . (23) 
-00 

For Gaussian stock fluctuations with the Hamiltonian ( 12T1) . the riskfree martingale mea- 
sure density has the form 16| 



p( M ™\x b ,t b ;x a ,t a ) =G(t b -t a 



-r w (t b -t a ) 



y/27Tv(t b - t a ) 



exp 



[x b X a Txw(tb ^a)] 
2v(t b - t a ) 



fx(t b )=x b f 1 pt b 

<3>(t b - * ) e -^(*»-*«) / Vx exp {-— [x - r Xw fdt \ . (24) 



x(ta ) — %a 



2v 



The corresponding time-compounded version of the superposition (TTTT) can be directly writ- 
ten as 



dvf^, t/ 5(v)P^ rw) (x b ,t b ;x a ,t a ) 



(25) 



Correctness of the latter can be verified by combining (|2"0"|) . (T2"2"|) . and (1241) . The result is [21 

poo 

P s (x b ,t b ;x a ,t a ) = / dv ft^,t/s(v) Pj; M ' rw) (x b ,t b ;x a ,t a ) 



-(r w t+Ax)/2(n M/6 



nT(t/S) 



\Ax — r w t\ 



1/2-t/s 



Ky2-t/s I Ax - r w t\ 



,(26) 
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with Ax = x b — x a and t b > t a . The measure density ( 126]) has for |Ax| ^> 1 the asymptotic 
behavior of the Erlang distribution 14| : 



exp 



\Ax\ 



Ax\ 



(27) 



which has a semi-fat tail (see Fig. [7]). It should be noted that the exponential suppression 
ensures that all momenta are finite. This fact is an important ingredient in showing that 
(126]) represents indeed the riskfree martingale measure density. The actual proof of the latter 
can be found in Ref. 16 . 




6 = 1,2, ji = 5 



t = 5 




FIG. 7: Normalized measure density (126j) : for various values of /i and 5 at fixed expiration time 
and (a 2 ) (left), for different expiration times at fixed 5 and fj, (right). We set the riskless interest 
rate r\y to be 12%. 



From the martingale measure (f26l) . we now calculate the option price at an arbitrary 
earlier time t a via the evolution equation 

/oo 
dx b 0(x b , t b )P s (x b , t b ; x a , t a ) . (28) 
-oo 

The value of the option at its expiration date t b is given by the difference between the 
underlying stock price S b = S(t b ) on expiration date and the strike price E, i.e., 

0(x b ,t b ) = Q(S b -E)(S b -E) = Q{x b -x E ){e Xb -e XE ) , (29) 

where xe = logE. Heaviside function Q(x) ensures that the owner of the option will exercise 
his right to buy the stock only if he profits, i.e., only when S b — E is positive. 

For Gaussian fluctuations of variance v , formula (1281) is evaluated with the measure ff24|) 
and yields directly the Black-Scholes option pricing formula [6]: 

Oi BS \x a ,t a ) = S(jt a )$(y+) - e-^-^E^y-), (30) 
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where represents the cumulative normal distribution and 



Vv 



\og[S(t a )/E] + (r w ± \v) (t b - t a ) 

\/v{t b - t a ) 



(31) 



The generalization to the present semi-heavy tail distribution is obtained from the superpo- 



sition 



23 



O s (x a ,t a ) = / dv ft„,t/s(v) 0[ BS \x a ,t a ) 



t > 0. 



(32) 



At this stage one should realize that integration in f[3"2"j) acts only on the &(y^) parts of 
Ov (x a ,t a ). This allows to write the option price in the Black-Scholes-like form, namely 



S (x a ,t a ) = S{t a )$M - e -rw(t b ~t a ) m (-) 



(33) 



where and ^ are Gamma-smeared versions of the functions and <&(y v ), re- 

spectively. For 5 — > 0, the new functions reduce, of course, to the un-smeared ones due to 
the limit ffT9]) . 

Let us now do the smearing operation for 



$(+) 



d^ W^W) = / dv f^ /6 (v) I -^t e^ 2 

</0 ./-oo V27T 



dv f tt i,t/s{v) 



di 



'0 ./-oo V27T 

By expressing the Heaviside function as 

dp 1 



-« 2 /2 



Q(x) 

and using Sokhotsky's formula 

1 



tpx 



p — i0_ 



P(-) + zvr5(p), 



(34) 



(35) 



(36) 



where V denotes the principal value of the integral, we perform the integral over £ and 
obtain 26 



$(+) 



1 

2 



V/2 



P 



(37) 



Here F^^ s (p, A) is the integral 



F t % s (p,A) = / dv/^e^ 
Jo 



(38) 
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This is evaluated as follows. We perform a change of variables from v to uj = Vvi, and 
define A = log[S(t a )/E] + r\yt, then we have 



,.t/S /-oo 

1 dcu to 21 / 5 - 1 e"^ 2 



nt/s) ,„ 

t/s 



exp 



%- \ uj H 

.2 V 



(39) 



In order to calculate gjP (p, A) we observe that the following integral can immediately be 
done: 



/•OO /*0O 

/ d/i// s -^W(p,A) = r(s) / dcu u 2t/5 - 2s - 1 exp 



.p ( 2A 

%- \ uj -\ 

2 V 



(40) 



Let us first assume that p > and Arg(A) = and perform the substitution to = y2A exp u. 
This yields directly 



2Acoshu-(2s-2t/<5)it 



poo POO 

/ duti'^gMfaA) = T(s)(2Af s - s du e» 

JO J-oo 

= l vT{s){e-™2A) t l 6 - s H^l 2t/s {pV2A) 



(p > 0) , (41) 



where Ha\z) is the Hankel function of the first kind 27]. This result holds for 3?(s) > 
t/5 — 1/2. Similarly we obtain for p < 0: 

POO 

/ d/i^-^ + )(p,A) = -i7rr( S )( e -2A)*/ 5 -^ 2t/a (-pV2A), (P<0)(42) 
Jo 



where ifa (z) is the Hankel function of the second kind [27J. The result is valid for 9ft(s) < 
t/5 + 1/2. 

Knowing these integrals, we find the function gj± \p, A) itself with the help of the Mellin 
inverse transform: 



c+ioo 



ds 
2^i 



ft 



i7iT(s)(e-™2Ay/ 5 - s H£l 2t/s (pV2A), p>0, 
-mT(s)(e™2A) t / & - s H { £_ 2t/5 (-pV2A), p<0 



(43) 



where c G (t/5 — 1/2, t/5 + 1/2). Inserting this back into (I3"9"|) and fl3T|) one can write for 
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$(+) 



t/«5 



dp e 



V/2 



r(t/5) y_ 00 27ri p 



(e~ i,r 2A)'/ 6 - s if, (1) 



2s-2t/<5 



(pV2A) 
-pV2A) 



1 //*/* 

2 + f(t75) 



C+lOO 



ds r(s)/x" s (2A) 



t/S-s 



dp e 



V/2 



o+ 2 ™ P 



L 2s-2t/5^ 



(44) 



Decomposing the Hankel functions into Bessel functions of first kind, 



this becomes 



$(+) 



#(1,2), 



(2A/x)*/ g 

r(t/(5) 



2 sin w 



(45) 



c+«oo 



ds r(s)(2A/x)" 



dp 



V/2 



27tz cos[7r(t/5 — s)] y p 

X J 2 s-2t/5(pV2A) + J2t/8~2s(pV2A) 



(46) 



The p-integral can now be easily performed yielding 



dp 



cos(ttC) i 0+ V I J-2c(pv2A) 

1 f (4) C r (l - i^i(C, 1 + 2C, -^4), 



ft c > o 



2^C 1 _ (a j c r (i + C ) i _ 2c, -A), » c < o 



(47) 



where ( = s — t/8. The function i-Fi(a, 6, is the Kummer confluent hypergeometric 
function 28|]. Note also that the fundamental strip c G (t/8 — 1/2, t/5 + 1/2) is for ( in the 
upper expression reduced to c G (0, t/<5+ 1/2), while for the lower expression c G (t/8— 1/2, 0). 
We now perform s-integral in (|4"6]) . and obtain 



with 

/4 +) (/i,t/5,A) 



= i + h^(fi,t/8,A) - h^(v,t/8,A) 



^rxc + t/5) (8/i) _ Cr(| _ c) iFi(C)1 + 2C _ A) 



(48) 



2T(t/6)y/t 



2r(t/5)0F 



Res 



E 

Res 



c 

r(c + 1/8) fA 2 fi 
C 



r(C + |) 1 F 1 (-C,l-2C,-A).(49) 
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To compute the residues of the poles of h\ and h 2 we need to decide in what way the 
poles are enclosed in the complex plane. Taking into account Stirling's large-argument limit 
of the Gamma functions: 



C 



t/S-3/2 



r(c + t/£)r(i + Q 
C 



2C 



( t/S ~ 



3/2 



(50) 



for \(\ — > oo (|Arg(£)| < tt), and the fact that the asymptotic behavior of iFi for large \(\ 
is given by Rummer's second formula [see, e.g., Ref. 29j Eq. (4.8.16)] 



1 F 1 (±(,1±2(,-A) « e~ A l\ |Arg(C)| < tt , 

iF 1 (C,l + 2C,-A) for C = -1,-2,-3,-4, ... undefined, (51) 



we obtain that the contour closure of h 2 depends entirely on the behavior of the Gamma 
functions. In fact, due to previous asymptotics we must close the contour in h 2 to the 
left as the value of the contour integral around the large arc is zero in the limit of infinite 
radius. Due to (8/i) _< > term hi closes the contour to the right. There are only simple poles 
contributing to h{ + \ which lie at ( = (2n + l)/2, n G N = (0, 1, 2, . . . ). If we assume for 
a moment that t/S ^ 1/2 + I, I G N, then the only singularities of h 2 are due to simple 
poles at ( = — 1/2 — m, m G N and ( = —t/S — k, k G N. Consequently we can write 



$(+) 



+ A 



+ 



2tt 



ra=0 




2_y>(t/5)_ n _ 1/2 (l/2) 



(l)2n+l V 2 



ra=0 



1 



2n+l 



l>5 



iFi ( n+-,2n + 2,-A 



cos(7rt/<5) \ Sv 



2A< 



2\t/S oo 



E 



(t/S), 



2A 2 \ n ft t 

J Q (l)2n+2t/«5(l)n \ Sv J 1 1 \ 6' 6' 



,(52) 



where we have set v = (v) — (a 2 ) and used the Pochhammer symbols (z)k = T(k + z)/T(z) . 
If 5 is very small, or t very large, we can use the asymptotic behavior (t/8) z — > (t/5) z and 
the fact that the sum in the third line of (1521) tends to zero due to strong suppression by 
l/(l)2n+t/<5- I n this case (1521 reduces to 



$(+) 



1 + \/z^ 1*1 ( 2' r 



1 3 {ytf 



Hvt), v = /3. 



(53) 



The calculation of is similar and is relegated to Appendix A. Here we state only the 
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result: 



2tt 



n=0 



4W _2_^(«/«)-n-l/2(V2) 




(1)2r+1 



n=0 



1 



2n+l 



iF x ( n + -,2n + 2,A 



+ 



cos(7rf/<5) \ <5?; 



2A 5 



2\t/S 00 



t>5 



2 \ n 



■ (54) 



In Appendix A we also show that for small 5, or large t, $i approaches the cumulative 
normal distribution $(y^ ). The asymptotic behaviors ensure us that the 5 + limit leads 
back to the original Black-Scholes formula, as it should. 

Another interesting situation arises for small A, where S(t a ) ~ Ee~ rw ^ th ~ ta ' . Using the 
fact that i-Fi(a, b, 0) = 1, and that we may neglect for small A the last two sums in Eq. (l52|) 
and Eq.( l54l) . we obtain that both and approach the cumulative normal distribu- 
tions $(y^))U=o and $(y^)U=o, respectively, implying that O s (x a ,t a ) -> (jB5) (x a , t a )U=o- 
Options with /I = are known as at-the-money-forward options, i.e. options whose strike 
price is equal to the current, prevailing price in the underlying forward market. Inasmuch, 
whenever the option is at-the-money-forward we regain back the formula of Black and Sc- 
holes. Since many real transactions in the over-the-counter markets are quoted and executed 
at or near at-the-money-forward 30J, this explains some of the empirical support for the 
Black-Scholes model . 

In this connection it is interesting to note that for put options, where the terminal con- 
dition is 



0{x b ,t b ) = Q(E - S(t b ))(E - S(t b )) 



(55) 



we would have obtained the pricing equation in the form 



O s (x a ,t a ) = Ee- rw ^\l-^)-S(t a )(l-^) 



(56) 



This shows that our option-pricing model fulfills important consistency condition known as 
the put- call parity relation [6j 



Of(x a ,t a ) = Of(x a ,t a )-S(t a ) + Ee- r ^-^ 



(57) 



where Of and Of denote put and call options, respectively. Consequently, the case when 
A « can be equally phased as a situation with Of(x a , t a ) ps Of(x a , t a ). We shall further see 
in the following section that the expansion factor A/\^j6 is directly related to the so-called 
moneyness of the options. 
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Moneyness is a measure of the degree to which an option is likely to have a nonzero value 
at the expiration date. In the Black-Scholes formula, the moneyness (measured in the units 
of standard deviation) is defined by 0] 

Vt + Vv \og[S(t a )/E] + r w t A 

mBS " ~ ^7t " TTt- (58) 

For at-the-money-forward options, the moneyness is zero. A positive value of m^s corre- 
sponds to in-the- money-forward options (i.e., options with positive monetary value) while a 
negative value represents out-of-the-money options. 

In our case v = a 2 is a random variable and for given t it is distributed according to the 
Gamma distribution ft^t/siv)- We may easily calculate the momenta of |m B s|- The odd 
moments are 



|m BS r +1 ) A 2n+1 (a" 2 "' 1 ) f A 2 \ n ftf (t/g)_ n _ 1/2 



(2n + l)! t«+V2 (2n+l)! \5v ) V Sv (l) 2n+1 



(59) 



so that the second sum in (152j) and (1541 corresponds to an expansion in the odd momenta of 
the moneyness. By taking into account that (t/6)- n = {t/5) n /{t/5 — n)2 n , the third sum in 
(|52|) and (I54p is seen to be an expansion in the even momenta of mes- Hence we can roughly 
characterize the three sums in the expansion (|52l) and ()54p as follows: the first sum is an 
expansion involving the properties of the stochastic process such as volatility, characteristic 
time, spread parameter 5 = q — 1), while the remaining two expansions are expansions 
involving the option contract characteristics (i.e., expiration time, strike price). 

The computations presented above assume a single-pole structure of (|46j) . i.e. that t/5 
is not half-integered. If t/5 is half-integered, the perturbation expansion (|52|) clearly fails. 
The half-integered t/5 case must be thus treated separately. For completeness, we state the 
corresponding result for $1 coming from the double poles in the s-plane of the integral 
( 146]) . First of all, the structure of stays the same because there are no multi-poles 
present. The corresponding computation of is not much more complicated than for 
simple-poles. Assuming that t/5 = 1/2 + 1, /gN, the analysis reveals that 



$(+) = - 
2 



v5 ^ (l/2 + /) n+1/2 (l/2) n / v5\ n _ / 1 




)2n+l 



iF 1 [n+ -, 2n + 2, -A 



/ A 2 ^ (l/2 + Q_ n _ 1/2 (l/2) B f 2A 2 Y „( 1 M0 A 



n=0 



I (-1)' / ^ ^ [V>(n-Z+1)+V>(n + 1)] /A 2 Y / 1 
+ r(/ + l/2)V2^^ r(n)r(n + /-l) iF 1 ^n+-,2n + 2 I -A 



,(60) 



where ^p(z) = T'(z)/T(z) is the Digamma function [31|. It is possible to check numerically 

that for t/5 — > 1/2 + /, / G N the relation fl52l) equals to (160]) (which is not true per- 

turbatively! ) . Thus is a smooth and non-singular function at aforementioned critical 
times. 
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V. CHARACTERISTIC TIME 



When confronted with a practical option price problem, one must decide whether or not t 
(i.e, the time to expiration) is large enough to deal satisfactorily with Gaussian distributions. 
This is done by estimating the characteristic (or crossover) time £*, below which the Black- 
Scholes formula is inapplicable and our solution becomes relevant. The estimate can be 
done with the help of a Chebyshev expansion j32|. We define a rescaled time-compounded 
variable 

Ax(t) — r T ,.,t 

z = y ' xw . (61) 

ytv 

Removing the trivial drift we find the difference between cumulative distribution at time t 



and that of the asymptotic Gaussian distribution as an expansion [32 



APx(u) 



dz 



P S (z,t) 



-u 2 /2 



'2rt 



-2 2 /2 



2tt 



Vi 



t 



Qiiu) Q 2 {u) 



+ p/ 2 



(62) 



Here we have used the abbreviation Pg(z, t) = P$(Ax — z, At — t), with the time t measured 
in the basic time units At - Functions Qj are Chebyshev-Hermite polynomials, the coeffi- 
cients of which depend only on the first j + 2 momenta of the random variable z appearing 
in the elementary distribution Ps(z,At ). The first two expansion functions are 



Qx{u) = 

QJ U ) = — u 5 + - — ) u + — -)u. (63 

v v 1 6! 8 V 3 9 J V 24 8 / 

Here k% and K4 are skewness and kurtosis of P$(z, At®). Since the drift was removed, 
Ps{z, Ato) is symmetric and the skewness k% vanishes. 

The characteristic time t* is now defined jfj] as a time at which the relative difference 
\APs(u)\/<&(—u) starts to be substantially smaller than 1 (to be specific we choose 1%), if 
u is taken to be a typical endpoint of the Gaussian central region, which we take as u = 1. 
Since \APs(u)\/<b(—u) has the value k 4 /(£24) for u = 1, we identify t* with k 4 . This implies 
that for t ^> t* = k 4 the Gaussian approximation is almost exact while for t < t* = k 4 the 
Black-Scholes analysis is unreliable and our new option pricing formula applies. 

To find k 4 we calculate the moment generating function G(p) associated with P$(z,At ) 
via the Fourier transform 



G(p) 



dy P(y, 1) 



(p 2 + 2 l i) l l s 



(64) 
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The cumulant generating function W(p) is then the logarithm of G(p). The cumulants c. 
are obtained from the derivativesand of W(p) : 



-0 



p=0 



S/i 



C2n 



(2n)\ 1 
n2 n 5fi n 



(2n) 
n2 n 



■if 1 6 



n-l 



C2n+1 



0, neN. 



(65) 



The excess kurtosis has the value k 4 = c 4 /(c 2 ) 2 = 35, implying that t* is equal to three times 
the spread of the distribution 8 — q — 1. This means that for large spread one can apply the 
Black-Scholes formula only a long time before expiration (long maturity options). 

It is worth noting that the present analysis is not applicable in cases when the distribution 
of returns is of the Tsallis type. This has heavy power-like tails, so that the second moment 
may be infinite and the above Chebyshev's expansion may not exist. In the present case 
the Tsallis distribution is in Fourier space and all cumulants are finite leading to the above 
estimate of t*. 

Our cumulant calculation (|65|) also further clarifies the meaning of the expansions (|52|) 
and (15^|) . In fact, the first sum in both (1521) and ([331) corresponds directly to the expansion 
in cumulants c 2n , n G N. This is because 



SC '"' ;i/2)„(f) • (66) 



Y(n)2 2n w ' \ 2 

Since v = j3, we can alternatively view the aforementioned expansion as an expansion in 
small P, 1.6. ? clS Sb ''high-Tsallis-temperature" -expansion. This establishes contact with the 
market temperature introduced in Ref. 
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Let us finally mention that the characteristic time can be crudely but fairly rapidly 
estimated from the shape of the variance (or volatility) distribution. In particular, the 
relative width of the distribution ft^t/siv) at time t is 



[v — v) 



5_ 
1 



(67) 



If the LHS of fl67|) is much smaller than 1 then the distribution is effectively 5-function 
and the volatility is a constant as assumed in the Black-Scholes analysis. However, when 
the relative width starts to be of order 1 the Black-Scholes model ceases to be valid. This 
happens at the time t* ~ 5 which is consistent with our previous estimate. 



VI. COMPARISON WITH EMPIRICAL MARKET DATA 



It is interesting to compare our solution (THHT) for European call options with realistic 
market data. Consider the option prices for an European option whose underlying is the 
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3/2/2003 3/2/2004 3/2/2005 3/2/2006 3/2/2007 

Time 

FIG. 8: Logarithmic plot of Dow Jones Euro Stoxx 50 index S(t) over 4 years (5 Feb 2003 - 21 
March 2007, in total 1057 trading days) with sampling interval Ato = 1 day. The index is fitted 
by a straight line, implying an exponential growth at an annual rate ~ 20%. 



Dow Jones Euro Stoxx 50 with the time series shown in Fig. [HJ The associated empirical 
option prices are plotted in Fig. [91 It is clear that because of the noise in the data no option 
pricing formula will fit the market prices perfectly. Even for a very good fit the option price 
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3/2/2003 3/2/2004 3/2/2005 3/2/2006 3/2/2007 

Time 

FIG. 9: Option prices for an European option with a strike price 4200 EUR and expiry on 21 
December 2007. The underlying of the option is the Dow Jones Euro Stoxx 50 with the time series 
given in Fig. 



would not be the most meaningful measure for the quality of the pricing model. Instead, one 
should test directly the hedging qualities of the pricing model j^] . To this end one introduces 
the A-hedge 
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Given A(r) we construct a so-called A-hedged portfolio n(r) by mixing stocks and options 
so that 



n(r) = O s (r) - A(r)5(r) . 



(69) 



The cash amount IT can be interpreted as the value of a portfolio of a trader who at time 
t bought one option Os(r) and sold an amount A(r) of the underlying with price S(t). In 
the ideal situation, i.e., neglecting the time delays in the determination of A(r) and the 
adaptation of the portfolio, and ignoring the transaction costs, this portfolio is perfectly 
hedged, i.e., it is free of fluctuations since the fluctuations of S and 0$ cancel each other |6j. 
The growth of the portfolio is therefore deterministic and proceeds at the riskfree rate r^: 



dn(r) 
dr 



r w U(r) 



(70) 



Any other growth rate would yield arbitrage possibilities. The amount A(r) deduced directly 
from empirical data should thus yield portfolio II (r) that is almost precisely e rwT (modulo 
multiplicative pre-factor). Consequently, A(r) computed from a good option pricing formula 
must give n(r) that is also close to the e rwT behavior. In our case the A-hedge is 



A = $l + \A) + 



d$i +) (A) 
OA 



S 



d^\A) 
OA 



Ee 



-r w t 



cL4 
dS 



d>«(A), (71) 



since the expression in the square bracket vanishes due to Eqs. (|48j) . (1B8I) and the identities 



[see e.g., Ref. [31 



dA 



~iFx(a + l,b + l,A) 
o 



S 



iFi(a,6,A) = e A 1 F l {b-a,b,-A) = j e rwt iF x (6 - a,b, -A) 
a iF 1 (a+ l,b,A) = a i-F\(a, 6, A) + A 



E 

dxF^a^A) 



OA 



(72) 
(73) 
(74) 



Let us note that by comparing (!69j) with the option-pricing formula (|33|) we can write 
the portfolio at the time r in the explicit form 



II(t) = -Ee- rw ^- r) &-\r) 



(75) 



Imperfectness of the hedging induced by the pricing formula ( 1531 thus depends on the 
actual behavior of $i (r) in time. In the Black-Scholes model the corresponding ^(y^) is 
effectively r-independent due to assumed form of the geometric Brownian motion for S(r). 
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In our model the situation is less obvious because the additional stochastic process due to 
volatility may substantially spoil the time independence. To see how big is the amount of 
the residual risk implied by our pricing model we turn now back to our empirical data. The 
best option-pricing solution O5 fit for the empirical option-price data from Fig. [9] is depicted 
in Fig. [TO! In the plot Fig. [10 (and the plots to follow) we have used Os calculated to 25th 
perturbation order with the double-poles removed. The obtained result is quite robust with 
the residual error smaller that 0.05%. 




3/2/2003 3/2/2004 3/2/2005 3/2/2006 3/2/2007 



Time 

FIG. 10: The best Os fit for the Dow Jones Euro Stoxx 50 option prices (cf. Fig. [9]). On this scale 
the Black-Scholes solution basically coincides with the Os prediction. 



Using as inputs: E = 4200 EUR, rw = 4.5%/year and tf, = 1270 trading days (t a = 
5 Feb 2003 = 0), the corresponding free parameters are then best fitted with 5 = 69.43, 
H = 351.29 (i.e., v = 0.000041). The fit can be further optimized when newly arrived option 
pricing data are taken into account. Details of the departure of the Os prediction from 
the Black-Scholes fit is depicted in Figs. [TTJ The goodness of the fit can be estimated, for 
instance, by the method of least squares. The correspondent likelihood function x 2 is f° r 
the Os fit in the period 21 Dec 006 — 21 Jan 2007 smaller by 1.2% in comparison with the 
Black-Scholes best fit. In the period 21 Jan 2007 - 21 Feb 2007 the difference in x 2 is 1.8%, 
and in the period 21 Feb 2007 — 21 Mar 2007 it is already 3.3%. This trend seems to get 
even more pronounced for periods closer to the expiry date. 

According to Section |V] the characteristic time corresponds to the time scale where non- 
Gaussian effects begin to smear out and beyond which the CLT begins to operate. For the 
call options at hand t* = 35 ~ 210 days, which, is particularly large and although the actual 
value will get further adjust with newly arrived option data, Fig. (TT] indicates that t* will not 
get dramatically changed. Note also, that despite the fact that t* = 210 days, the departure 
effect is visible already 10 months before the maturity. 

By having the best Os fit we can construct the daily portfolio II according the prescription 
( |69|) . The resulting portfolio is shown in Fig. [121 To quantify the fluctuations it is convenient 
to write U(r) = -E^i~\r) e-^ b ~^ rw = -E$ (y~) e -^ b ~ T ^ rw+Srw ^ . The relative interest- 
rate fluctuations 6r w /r w implied by the stochastic nature of the volatility are according to 
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2/21/2007 3/7/2007 3/21/2007 

Time 

FIG. 11: Two successive details of the departure of the 0$ fit from the Black-Scholes values. Black- 
Scholes option prices are taken with v = v. The roughly 5% departure is clearly visible already 
9 — 10 months before the expiry date. 



Fig. [12j 0.007(±0.002). Consequently we can conclude that for the data at hand the 0$ 
pricing formula induces a A-hedging strategy that is close to being optimal. 

VII. CONCLUSIONS AND OUTLOOK 

We have developed the theory of an option pricing model with a stochastic volatility fol- 
lowing a Chi distribution. The corresponding volatility variance is then distributed with an 
ubiquitous Gamma distribution. Our direct motivation was drawn from the high-frequency 
S&P 500 volatility fluctuation data, whose distribution is well fitted in this way. There are 
two interesting implications resulting from such volatility fluctuations. First, the returns 
of the corresponding asset prices exhibit semi-heavy tails around the peak (i.e., leptocurtic 
behavior) which mimics empirically observed long-range correlations. At the same time, 
the returns preserve some of typical features of the original Black-Scholes model, namely 
they follow a linear ltd stochastic equation with multiplicative noise (though with stochastic 
volatility) and a continuous stock dynamics. Second, the associated density operator in 
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FIG. 12: Negative daily portfolio —II constructed from the O5 pricing formula (|33p through A- 
hedging prescription (I68p and (I69h , The best exponential fit confirms the annual interest rate 4.5%. 
The relative fluctuations br w {r)jr w = 0.007 ± 0.002. 



momentum space is of the Tsallis type. The inverse Tsallis temperature then agrees with 
the average variance of the distribution. 

Our main result is a generalized Black-Scholes pricing formula that takes into account 
the above volatility behavior. With the help of a Mellin transform we were able to find 
an ensuing analytic solution for the price of European call options. The result is expressed 
as a series in the higher normalized cumulants and in higher moments of moneyness. Due 
to a spread parameter 6 related to the extra kurtosis of the log-return data, our model is 
capable of calibration to a richer set of observed market histories than the simple Black- 
Scholes model, which is a special case corresponding to a zero spread 5, or to a very long 
time horizon (t ^> t* ~ $). Comparisons with other time-dependent volatility models such 
as ARCH [i], GARCH 0j or multiscale GARCH [4] will be addressed in future work. 

In the light of recent works on superstatistics we should mention that the density operator 
representation (fT5l) reveals the sup erstatistic nature of the Tsallis distribution in momentum 
space. Recently, C. Beck [33|, |34j and C. Beck and E. Cohen [35j, prompted by the works 
by G. Wilk and Z. Wlodarczyk 36(, have suggested that the origin of certain heavy-tail 
distributions should be understood as weighted averages of the usual exponential statis- 
tics. Such averages have been used also in the textbook [16| to calculate option prices for 
non-Gaussian distributions 
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Recently, F. Sattin [37J rephrased the same procedure in 
terms of evolving systems embedded within a static but non-trivial background. All these 
approaches commonly strive to interpret broad distributions as a result of an averaging of 
usual exponential (Gibbs) distributions over certain fluctua ting (random) parameter. This 
is the procedure running under the name superstatistics [35j, |37| • Although frequently used, 
superstatistics procedure is still far from being systematized and, in fact, smearing distri- 
butions are usually chosen to fit experimental data. The theory of option pricing in this 
work may therefore be viewed as an example for superstatistics with a smeared Gaussian 
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representing the probability distribution of volatilities. 

Finally, the present scenario also provides an interesting meaning to the Tsallis escort 
parameter q. In option pricing models where stock fluctuations are directly fitted by a 



Tsallis distribution of returns 16|, |l8j], the q parameter is a fit parameter with a typical 
value around 1.5. In the present case where the momentum distribution is of the Tsallis 
type, q — 1 is proportional to the characteristic time t* below which the Gaussian treatment 
of stock fluctuations is inadequate. 
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APPENDIX A: FOKKER-PLANCK EQUATION FOR STOCHASTIC EQ.(l) 
AND SOME CONSEQUENCES 

Some remarks may be useful concerning the stochastic equations for volatility and vari- 
ance stated in the introduction. 

Let us assume that the variance v is driven by the Ito stochastic process 

dv(t) = i{t)[v(t) -fi(t)v(t) - a(v(t),fi(t),u(t))]dt + ^2^{t)v{t) dW{t) . (Al) 

Here j(t),fj,(t) and u(t) are arbitrary non-singular positive real functions on IR + . The func- 
tion a(v(t), (j,(t), v(t)) is assumed to be a non-singular function of its arguments. The corre- 
sponding Fokker-Planck equation for the distribution function p(v, t) reads 

9 -P^ = ^-Mt)[v^t)-u(t) + a(v^(t),u(t))]pM] + ^h(t)v P (v,t)]. (A2) 

Let us choose the function a(v,ji(t), u(t)) to be 
(log/i)' 



7 



a(v,n,u) = v^p- + ue ""{[r(i/)log(^) +r(z W )(^z/)-logM)] {vnY v 

(A3) 



v 2 J 7 



where T(x,y) is the incomplete gamma function, ip(x) the digamma function, and 2-^2 a 



hypergeometric function 3JJ) the solution of (IA2I) has the form 



p(v, t) = f ^ ^VH-^ . (A4) 
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This is precisely the Gamma distribution f^(t)^(t)iy) in Eq. <^ describing the temporal 
statistical evolution of the random variable v. Note also that 7(t) does not enter in the 
solution. 



By applying the Ito formula of the stochastic calculus 39( one easily obtains the corre- 
sponding Ito stochastic equation for the volatility o = ^fv in the form 



dcr(t) 



7(*) 



(y(t) - 1/2) _ _ a(a 2 (t),p(t), u(t)) 



a(t) 



a{t) 



df - \/2&ldW(t), (A5) 



from which we derive the Fokker-Planck equation for the distribution of the volatility a 



dp(a, t) 
dt 



+ 



<9a \ 2 

7(*) 



a 2 
9^ 



4 



(Tp(t) 

p(a, t) 



K*)-l/2) , a((7 2 , //(*),!/(*)) 



(7 



(7 



p((T,t) 



(A6) 



Of course, this equation is simply related to ( 1A2I) by the substitution v = a 2 . Equation 
Eq.( 1A6l) has the solution 



p(a,t) = 2af mAt) {a 2 ) , (A7) 

which is also known as the Chi distribution 

The particular solutions ( 1A4I) and 1 ) A 71) are the only ones fulfilling the asymptotic con- 
ditions p(v,t)\ t ->oc = 8(v — v) and p(cx, i)|t-»oc = S(a — a), respectively. They are also both 
normalized to unity. Note that the uniqueness of the solutions ( 1A4I) and ( 1A7I) ensures that 
the stochastic equations ( 1A1I) and ( 1A6I) have weak solutions. 

Anticipating the empirical form presented in Section [III we shall now assume that at large 
t the functions p[t) and u(t) behave like pit) ~ pt and z/(t) « z/t. Without loss of generality 
we also set j(t) = 1. In this case the large-t behavior of a(v,p(t), u(t)) can be easily found. 
Using the asymptotic formulas [29( 

^(^)k-oo ~ log(i/t), 

2F 2 (z/t,z/t;z/t + l,z/t + l;-^)| t ^ 00 « e^* (l + , 

r(z/t,t; / ut)|^ 00 « e^i/*)"*- 1 / 2 ^ - (i^M) 1 ^-^' (l + ^) , (A8) 



we obtain that a(v,p(t), v(t))\t-+oo ~ 0(l/t). In this limit, the stochastic differential equa- 
tion flAlj) for the variance corresponds to the Cox-Ingersoll-Ross process |38(. Although this 
asymptotic equation resembles Heston's stochastic volatility model [lj], there is a difference: 
in the drift term v(t) and p(t) are linear functions of t rather than constants. 
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Let us assume the same linear behavior fi(t) m fit and v(t) ~ id at small t (more precisely 
for id >C 1). For convenience we set fi/u = v, as in Section HTT1 The small-t behavior of 
a(. . .) can be easily found from the asymptotic formulas [29| 



2F 2 (ist,irt;vt + l,irt + l;—v/vvt)\ u t-to ~ 1, 

r(i4,v/v 1^)1^0 ~ -log[i/t]-log[v/7j]-C, (A9) 

where C = 0.57731... is the Euler-Mascheroni constant. Then a(v, fi(t), v(t))\vt->o ~ 
—v/(vtv) + 0(1/ vt). This shows that in the small- 1 limit, the drift term is dominated 
by the function a(v,fi(t), v{t)). Let us finally note that the stochastic equation (IA5I) corre- 
sponds to an additive process, rather than to a multiplicative one. 



APPENDIX B: CALCULATION OF $ 



(-) 



Let us calculate the function $1 ^in Eq. (I54p appearing in the generalized Black-Scholes 
formula (1331) . The procedure is analogous to that for We start from Eq. (1371) which 

reads now 



$(-) 



1 ^r^w f 



-p 2 /2 



27TZ **»>*/* ^ p 



where 



t/5 roc 



*C/,(p) = 2 



r(t/5) 



/•oo 

/ do; to 2 '/ 5 - 1 e"^ 2 exp 
Jo 



p (2A 

i- uj 

2 V uj 



= 2 



t/s 



T(t/6) 



' p (2A 

i— uj 

2 I w 



-ipv^Acoshu— (2s— 2t/(5)« 



(Bl) 



(B2) 



To find ^(p, A) we first calculate the integral 

poo poo 

/ d /U/ x s - 1 ^-)(p,A) = T(s) do; a; 2 */ 5 - 2 - 1 exp 
Jo Jo 

/oo 
dw e" 
-oo 

= 2r( S )(e- i7r 2A)'/ <5 - s ^ 2t/(5 (pv / 2l), (p>0), (B3) 



where ifa(z) is the modified Bessel function [27]] . The result (1B3P holds in a strip t/5 — 1/2 < 
3?(s) < t/<5 + 1/2. Similar calculations for p < yields 

/>oo 

/ dfifi-^faA) = 2T{s){e™2A) t / 5 - s K 2s _ 2t/s {-pV2A) 1 (p < 0) . (B4) 
Jo 
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Again, the result is true in a strip t/8 — 1/2 < 3?(s) < t/5+ 1/2. With the help of the Mellin 
inverse transform we now find 



- )(M)= r^^4 ,rtf ' , "' ww ' p> °-^ <B5) 

' 2m p 'I (rtif -Jfj.^.l-pv^), p<0, ' 



C— lOO 



where c G (t/5 - 1/2, t/8 + 1/2). Inserting this back into (|B2|) and (jBljl we obtain 

v ~2 + T(t/8) J _ tori p J ,_j 0O 7rz ^ (s) j ( e - 2 A)^-^ 2s _ 2t/5 (-pV2A) 

= ? + t^ttft / 5- r( S )(2A/,)^- s sin[7r( S - t/5)] / -2— -K 2s _ 2t/J (pV2l) . 
2 r(t/d) y c _ ioo 2vn J 0+ tt p 

(B6) 

The p- integration can be carried out explicitly yielding (£ = s — t/8) 



sin(TrC) / dp K 2( (pV2A) 

J0-L. V 



sin(7rC) 



'7T 

8^ 



[A- f r(-c)r(20 iF!(-c, i - 2C, A) + ^r(c)r(-2C) 1*1 (C, 1 + 2c, a 
( -J r(i - c) iiMC, i + 2C, A) + (~) ? r(c + ±) M-c, 1 - 2c, A) 



• (B7) 

The integration is valid for 3? £ > and 9ft £ < 0. The result allows us to write as 

$H = ~ - h[-\^t/8,A) - h^(ji,t/5,A) , (B8) 

with 



fcTVt/M) = 2r(f/f)A E ( ' (8ri- g r(i-c) 1 Fi(C.i + 2C,A), 
*'<"^> = S^E^(^)^r(C + i)^W,i-2CM). 

(B9) 

Mellin's fundamental strip for £ can be conveniently chosen in as < £ < 1/2, and 
in h 2 ~^ as — 1/2 < £ < 0. As previously in the calculation of the contour of the 
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s- integration for ls closed on the left. Thus we obtain by analogy with $1, that the 
final expression has the form (1541) . 

In the small-5 limit where (t/8) z — > (t/5) z , the third sum in (1541) goes to zero. In such a 
situation 

$(-) = - 
2 

i.e., we regain the cumulative normal distribution. In the limit A — > 0, one finds that 



1 + 



2/s i^i 




(BIO) 



APPENDIX C: INFORMATION ENTROPY AND TSALLIS STATISTICS 

A useful conceptual frame that allows to generate important classes of distributions is 
based on information entropies. Information entropies generally represent measures of un- 
certainty inherent in a distribution describing a given statistical or information-theoretical 
system. Central role of information entropies is in that they serve as inference functionals 
whose extremalization subject to certain constraint conditions (known as prior information), 
yields the MaxEnt distribution. Importance of information entropies as tools for inductive 
inference (i.e., inference where new information is given in terms of expected values) was 



emphasized by many authors 40 



Among the many possible information entropies one may focus attention on two examples, 
first on Renyi's entropy 

sf) = -^io g J>?, g>o, (Bii) 

i 

and second on the Tsallis-Havrda-Charvat (THC) entropy 

S f HC) = J^- q fe^" 1 ) ' g> °- (B12) 

The discrete distribution V = {pi} is usually associated with a discrete set of micro-states 
in statistical physics, or set of all transmittable messages in information theory. In the limit 
q — > 1, the two entropies coincide with each other, both reducing to the thermodynamic 
Shannon-Gibbs entropy 

S = -J^PilogPi. (B13) 

i 

Thus the parameter q — 1 characterizes the departure from the usual Boltzmann-Gibbs 
statistics or from Shannonian information theory. 

It is well known [4]J that within the context of Shannonian information theory the laws 
of equilibrium statistical mechanics can be viewed as inferences based entirely on prior 



30 



information that is given in terms of expected values of energy and number of particles, 
energy and volume, energy and angular momentum, etc. For the sake of simplicity we 
shall consider here only the analog of canonical ensembles, where the prior information is 
characterized by a fixed energy expectation value. The corresponding MaxEnt distributions 
for 5 W and sf HC) can be obtained by extremizing the associated inference functionals 



L iTHC) {v) 



-^—log^pl - a^pi - P(H) r , 

i i 



(B14) 



where a and j3 are Lagrange multipliers, the latter being the analog of the inverse temper- 
ature in natural units. The subscript r on the energy expectation value (H) distinguishes 
two conceptually different approaches. In information theory one typically uses the linear 
mean, i.e., 

(H) 1 = (H) r=1 = Y,P* E *> (B15) 

i 

while in non-extensive thermostatistics it is customary to utilize a non-linear mean 



(H) q = (H) r=q = with P^q) 



Pi 



q • 



Pi 



1 . 



(B16) 



The distribution Pi(q) is called escort or zooming distribution and it has its ori gin in chaotic 
dynamics 42j and in the physics of multifractals [43]. Simple analysis reveals 44j that 



5Lf\y) 



P 



(i) 



1/(9-1) 



P 



(2) 



1 - P(q - l)AEi 
Z£[l - P(l-q)AE i } 1 ^ 



for (H) r=1 , 
for (H) r=g . 



(B17) 



Here f3 = f3/q and AEi = E^ — (H) r . By the same token one obtains for the THC case 44 



5Lf HC) 



(V) 



P 



(i) 



7 -i 

THC 



1/(9-1) 



o~Pi 



P 



(2) 



1 - p*(q-l)AEi 
Z^ HC [l - P*(l-q)AE l ] 1/{1 - q) 



for (H) 
for (H) 



r=l 



(B18) 



r=q i 



with (3* = (3/J2iPi and $* = P/J2iP 9 i- So in contrast to (JBT7]), the THC MaxEnt 
distributions are self-referential. Generalized distributions of the form (1B17I) and (1B18j) are 
known as Tsallis distributions and they appear in numerous statistical systems [451 ] . For 
historical reasons is = {p[^} in (1B18[) also known as the Bashkirov's 1-st version of 
thermostatistics, while = {p^} in ( 1B18I) is called the Tsallis' 3rd version of thermo- 
statistics. An important feature of Tsallis distributions is that they are invariant under 
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uniform shifts e of the energy spectrum. So one can always choose to work directly with E{ 
rather than AE%. 
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